3 Décembre 2024

Plan

  • Partie 1: Analyse du transcriptome par RNA-seq

    • Pipeline RNA-seq

    • Analyse en Composantes Principales (ACP)

    • Analyse Fonctionnelle, enrichissement GO

    • Quelques banques publiques utiles (Génomes, annotations, ontologies)

  • Partie 2: scRNA-seq

    • Pipeline

    • Analyse différentielle

    • Intégration de plusieurs échantillons

Partie 1: Analyse d’expression des gènes: RNA-seq

Je me présente

Qu’est ce que la bioinformatique ?

  • Apparition en 1970: B Hesper et P Hogeweg, « Bioinformatica: een werkconcept », Kameleon, vol. 1, no 6, 1970, p. 28–29

La bio-informatique est constituée par l’ensemble des concepts et des techniques nécessaires à l’interprétation informatique de l’information biologique. Plusieurs champs d’application ou sous-disciplines de la bio-informatique se sont constitués (Wikipedia):

  • La bio-informatique des séquences
  • La bio-informatique structurale
  • La bio informatique des réseaux
  • La bio-informatique statistique et des populations

Projet du séquençage du Génome Humain

  • Idée lancée en 1985 par 3 scientifiques, Renatto Dulbecco, Robert Sinsheimer (directeur de UCSC) et Charles DeLisi, qui financera le projet (Directeur dept. de biologie du Département de l’Energie US)
  • Séquençage lancé en 1988 par le National Research Council. En suisse est créé HUGO (Human Genome Organisation) pour la coordination.
  • En 1998, Craig Venter, crée Celera Genomics avec pour objectif le séquençage en 3 ans par séquençage Shotgun et le brevetage du génome (!).
  • En 2000, la complétion du séquençage est annoncé pour le Consortium et Celera Genomics (match nul) par le président B. Clinton. Coût: $3B.
  • Celera Genomics avouera avoir utilisé les données du Consortium pour son propre assemblage, mais reproduira un séquençage de Novo 3 ans plus tard…
  • Publication des séquences brutes en 2001 et des séquences finales en 2004.

Différentes technologies de séquençage

Déroulé d’un séquençage

Le principe d’un séquençage NGS consiste à:

  • Créer une librairie de fragments d’ADN (par fragmentation enzymatique et mécanique de l’ADN génomique).
  • Relier ces fragments à des adaptateurs, des petites molécules d’ADN de séquences connues nécessaires au séquenceur.
  • Une sélection de la taille minimum et maximum des fragments est effectuée pour des raisons techniques, notamment l’amplification PCR.
  • Faire une amplification PCR des fragments.
  • Faire le séquençage.

Principe général du séquençage

Séquençage Shotgun

Génération des séquences (Illumina)

Séquençage par Illumina - principe

  • Hybridation d’un brin sur un oligonucléotide attaché à la FlowCell

  • Un brin complémentaire est synthétisé.

  • La molécule d’origine est enlevée et la molécule libre s’hybride en pont.

  • Un brin supplémentaire est synthétisé de nouveau.

  • Les brins complémentaires à la cellule d’origine sont lavés et il ne reste que plusieurs copies d’une même brun (clusters)

  • Il reste à séquencer les brins présents: Lors de cette étape, le nucléotide incorporé est identifié grâce à un groupe fluorescent identifié par laser, permettant d’enregistrer la séquence de manière informatique.

  • https://www.youtube.com/watch?v=CZeN-IgjYCo

  • https://www.youtube.com/watch?v=WneZp3fSJIk

Figure B: Schéma représentant les clusters sur une FlowCell

Figure C: Réaction de séquençage

Enrichissement

Il est possible de créer une librairie enrichie en régions d’intérêt, par exemple pour séquencer uniquement les régions codantes du génome:

  • Par capture: Il est possible de concevoir des sondes qui vont s’attacher à l’ADN d’intérêt, elle même étant liées à des molécules de biotines attachées à des billes magnétiques, qui sont conservées après lavage.
  • par amplicon: Il est possible de faire une amplification sélectionnée par des amorces PCR choisies.

Séquençage par capture

Séquençage par amplicon

Différence Pair-End (PE) et Single-End (SE)

Il est avantageux d’obtenir des fragments de lecture les plus longs possible pour un alignement le plus fiable possible.

  • Lors d’un séquençage simple (Single-end), les brins sont séquencés en partant d’un unique adaptateur. Un read correspond donc ensuite à un fragment.

  • Lors d’un séquençage appairé (Paired-end), les brins sons séquencés à partir de leur deux extrémités. Les fragment résultats sont appelés R1 et R2 et sont liés, qu’ils soient recouvrant ou pas.

Un site expliquant cela de manière intéressante:

http://thegenomefactory.blogspot.com/2013/08/paired-end-read-confusion-library.html

Pair-End (PE) et Single-End (SE)

Analyse RNA-seq - principe

Analyse de l’expression des gènes = le Transcriptome. C’est une grandeur dynamique.

  • Technologie à haut débit précédente: les Puces à ADN.
  • Technologie basée sur le NGS: Le RNA-seq.

Beaucoup d’outils sont communs à ces technologies.

Le NGS appliqué à l’analyse du transcriptome permet:

  • Une meilleure concordance entre plateformes
  • Une forte sensibilité et meilleure dynamique
  • Toutes espèces, toutes régions transcrites
  • Une variété d’applications (Analyse diférentielle, analyse de l’épissage alternatif, analyse des gènes de fusion, séquençage de novo)

Mais…: Complexité et coût calculatoire accrus = pipeline bioinformatique plus complexe par rapport aux microarrays.

Déroulement d’une analyse RNA-seq

  • Séquençage
  • Contrôle qualité
  • Alignement sur génome de référence
  • Quantification des valeurs d’expression (comptage) et normalisation
  • Analyse différentielle entre conditions expérimentales
  • Visualisation des données (“diagrammes de chaleur”, “volcano plots”).
  • Analyse fonctionnelle des gènes (“Enrichment Analysis”)

L’analyse bioinformatique fait partie intégrante du processus de traitement.

Analyses possibles par le RNA-seq:

  • Assemblage de transcriptome de novo
  • Analyse de l’épissage alternatif
  • Découverte des gènes de fusion
  • Découverte de nouvelles classes (exemple: sous type moléculaire de tumeurs) (analyse non supervisée)
  • Analyse différentielle: C’est cette application que nous allons décrire ici.

Principe du RNA-seq (1)


Source: Thomas Girke, Analysis of RNA-Seq Data with R/Bioconductor, December 14, 2013)

Principe du RNA-seq (2)


Source: Wikipedia

Départ: les fichiers issus du séquenceur (Fichiers FASTQ)

Ils contiennent les reads: petite séquence d’un fragment d’ADN de longueurs plus ou mons fixe.

  • Single-end
    • Chaque read est indépendant
  • Paired-end
    • Le séquençage est fait par chaque extrémité de chaque brin. Dans ce cas, les reads sont organisés par paires.
@HWI-ST865:166:D0C4KACXX:2:1101:1042:1954 1:Y:0:
CNANAAATNAANNNNGNNNNNNNNNANNNNNAAANNNTNNNNNNNNNTNNTGNNNNTTGTTTNNTTGTGGGTTTCTCTGTCCCCN
+
#####################################################################################
@HWI-ST865:166:D0C4KACXX:2:1101:1241:1970 1:N:0:
CCAGCGACACTTGCAGCTTAGGGGCAAGAGGCTCCCACAACACCCTGTGCGATCGGAAGAGCGGTTCAGCAGGGATGCCGCGGCC
+
GFFIGIIIFGEHHIJJJIIGGGHIIBD=BFG?EDECC@FGCHC?BCCBB)53(;;B;?8299?######################

Mesure et encodage qualité: le Phred

Parenthèse: Infrastructure pour le stockage et le traitement des données

Il est important de stocker les données dans un emplacement fiable, avec beaucoup de mémoire vive et d’espace disponible, et accessible par un réseau rapide (Bref, un serveur de stockage).

Pour le traitement des données, un cluster de calcul sous linux est idéal, une machine Linux ou MacOS convient parfaitement.

Il existe également des logiciels commerciaux sous Windows.

Contrôle Qualité par FastQC

Alignement par STAR: Principe

STAR est un aligneur extrêmement précis qui peut dépasser certains aligneurs par un facteur 50 de vitesse. Par contre, il est très gourmand en mémoire.

Son principe de fonctionnement est basé sur deux étapes principales:

  • Recherche de graines (seeds)
  • Groupement (clustering), agrafage et calcul de score

Source: https://hbctraining.github.io/Intro-to-rnaseq-hpc-O2/lessons/03_alignment.html

Principe de STAR

Etape 1: Recherche d’alignement exact maximum pour chaque read (Maximal Mappable Prefixes - MMP) pour obtenir une liste de graines (seed1)

Principe de STAR

Etape 2: STAR va ensuite rechercher un alignement exact pour les portions non alignées des graines (appelées seed2).

Principe de STAR

Etape 3: Si l’alignement exact n’est pas possible, STAR va étendre la séquence recherchée.

Principe de STAR

Etape 4: Si l’alignement n’est toujours pas possible, STAR va couper les séquences adaptatrices ou contaminantes ou de faible qualité.

Principe de STAR

Etape 5: Groupement, Agrafage et Calcul de score. Les graines séparées sont groupées et agrafées sur la base sur leur proximité avec un site donneur pour créer un read complet.

Visualisation par IGV (Integrated Genomics Viewer)

IGV est un explorateur de génome. Il permet de visualiser de l’information génomique le long du ou des chromosomes de l’organisme que nous sommes en train d’étudier.

  • La première étape est de choisir un génome de référence (Ex: hg19 ou mm10).
  • Une boîte permet de situer la position actuellement visualisée sur le chromosome.
  • Les informations correspondant à chaque échantillon sont présentées sous forme de pistes (tracks).
  • Les gènes sont affichés dans leur propre piste (annotation) avec les exons et introns. Il est possible d’afficher les transcrits individuellement.
  • Plusieurs types de signaux peuvent être affichés (fragments de lecture, segments CGH, couverture, etc…).

Exemple d’alignement (Mapping) RNA-seq sur génome de référence

Analyse de séquence avec IGV

Quelques définitions: Couverture et profondeur

Annotation du génome de référence

Les annotations du génome de référence sont disponibles sous forme de fichiers GFF/GTFou BED auprès de Ensembl (BioMart - https://www.ensembl.org/info/data/ftp/index.html) ou NCBI (https://www.ncbi.nlm.nih.gov/refseq/).

#!genome-build GRCh38.p13
#!genome-version GRCh38
#!genome-date 2013-12
#!genome-build-accession NCBI:GCA_000001405.28
#!genebuild-last-updated 2019-06
1   havana  gene    11869   14409   .   +   .   gene_id "ENSG00000223972"; gene_version "5"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene";
1   havana  transcript  11869   14409   .   +   .   gene_id "ENSG00000223972"; gene_version "5"; transcript_id "ENST00000456328"; transcript_version "2"; gene_name "DDX11L1"; gene_source "havana"; gene_biotype "transcribed_unprocessed_pseudogene"; transcript_name "DDX11L1-202"; transcript_source "havana"; transcript_biotype "lncRNA"; tag "basic"; transcript_support_level "1";

Comptage

  • Sous l’hypothèse que le nombre de reads venant d’un certain gène est proportionnel à l’abondance de son ARN dans la cellule, on compte les reads venant de chaque gène, transcrit ou exon du génome.

  • Il est possible de faire un script ‘maison’ mais il existe maintenant un grand nombre de programmes pour faire cette fonction, notamment featureCount de la suite logicielle SubRead.

  • Les localisations génomiques des objects génomiques (genome features) sont données en entrée du programme de comptage (Fichier d’annotation GTF/GFF), permettant d’assigner les comptages à chaque transcrit, gene, ou exon.

Comptage par featureCounts

featureCounts permet de compter les fragments de lecture correspondant à chaque gène et d’obtenir un tableau récapitulatif pour l’ensemble de nos échantillons.

featureCounts du package subRead permet de faire ce comptage à partir des fichiers BAM résultants de l’alignement de séquences par STAR et un fichier d’annotation de génome au format GFF.

Il est possible de travailler en mode PE (Paired-End) ou SE (Single-End). Il n’est pas obligatoire de trier les BAM avant de les donner à featureCounts.

Assignation des reads

Les fragments de lecture sont assignés à des attributs (feature) ou des meta-attributs (meta-features). Il y a donc 2 modes.

  • Les attributs (features) correspondent à une zone contiguë de l’ADN, par exemple un exon.

  • Les meta-attributs (meta-features) correspondent à un objet biologique résultant de l’assemblage de plusieurs attributs (les gènes)

featureCounts peut compter au niveau des exons ou des gènes. Il est recommandé d’utiliser des identifiants uniques pour les meta-attributs, comme par exemple les geneID du NCBI, et d’éviter les noms de gènes.

Alignements multiples et recouvrants

Pour les alignements multiples, featureCounts peut les traiter de 3 manières:

  • Les ignorer (Défaut pour la version bash)
  • Les compter plusieurs fois
  • Les compter plusieurs fois mais avec une contribution calculée en fraction (\(1/n\), \(n\) étant le nombre de fois ou le read apparaît).

Pour les alignements recouvrants plusieurs attributs, featureCounts peut les traiter de 2 manières:

  • Les ignorer (Recommandé pour le RNA-seq)
  • Les compter pour chaque attribut

Obtention du comptage pour les transcrits

featureCounts permet de compter les reads au niveau des exons (features) ou des gènes (meta-features). Par contre, Il n’est pas possible d’obtenir les reads pour les transcrits, car il ne peut déterminer à quel transcrit les assigner.

Ce type d’analyse est considérablement plus complexe. Pour ce faire, il faut utiliser un pipeline spécifique basé sur l’aligneur Salmon, qui permet d’assigner les reads aux transcrits sur une base de probabilité.

Voir la documentation pour l’aligneur Salmon: https://salmon.readthedocs.io/en/latest/salmon.html.

Normalisation RPKM

Nous avons différentes tailles de libraires, des biais de séquences dus à la PCR, ou un contenu RNA qui peut différer en échantillons. Il est donc nécessaire d’appliquer une normalisation pour rendre les échantillons comparables.

Méthode la plus utilisée: RPKM: Reads Per Kilobase per Million. Soir une valeur de comptage \(read(G)\) pour un gène:

  • On divise la taille de la librairie par \(10^6\). Cela nous donne le facteur par million.
  • On divise la valeur de comptage \(read(G)\) par ce facteur pour obtenir des reads per millions (RPM).
  • On divise la valeur RPM par la longueur du gène en Kilobases, pour obtenir des RPKM.

\(RPKM(G) = \frac{[\sum read(G)]}{[\sum Read].longueur(G)}.{10}^{6}.{10}^{3}\)

Variante: FPKM pour le Paired-end. Dans ce cas, deux reads peuvent appartenir au même fragment et ne sont pas comptés 2 fois.

Effet des gènes très ou peu exprimés sur la normalisation RPKM

Il faut voir une expérience RNA-seq comme un échantillonnage de l’espace des transcrits par la librairie.

Une expérimentation RNA-seq donne donc la proportion des transcrits pour une taille de librairie donnée.

Si un gène est fortement exprimé dans un échantillon et peu exprimé dans l’autre, on aura une augmentation artificielle de l’expression de l’ensemble des gènes (!).

Gène Longueur Count Ech1 Count Ech2 RPKM Ech1 RPKM Ech2
Taille librairie 20 20
Gene 1 2 8 0 0.2 0
Gene 2 2 6 10 0.15 0.25
Gene 3 2 6 10 0.15 0.25

Normalisation TMM

TMM - Trimmed mean of M-values. C’est l’une des normalisations recommandées qui permet d’éviter l’effet des gènes fortement ou faiblement exprimés (Robinson & Oshlack. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology volume 11, Article number: R25, 2010).

Au lieu de faire une mise à l’échelle propre à une librairie, un facteur de normalisation global est calculé en assumant que la majorité des gènes n’est pas différentiellement exprimée et en ne tenant pas compte des valeurs extrêmes de comptages.

Ce facteur est calculé à partir de l’abondance relative de l’expression de chaque échantillon, soit un fold change global (Ce que l’on appelle les M-Values).

Un échantillon est choisi comme référence. L’abondance relative moyenne est calculée en “trimmant” les distributions de logFC (d’où le Trimmmed Mean), ce qui donne un facteur correctif pour chaque échantillon.

Normalisation TMM

Il est calculable par le package d’analyse edgeR sous R/Bioconductor.

Autre normalisation: celle proposée par DESeq2.

Dillies et al. A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Brief Bioinform . 2013 Nov;14(6):671-83.

Normalisation (Suite)

La Normalisation type TMM ou autre n’est nécessaire que pour la visualisation ou l’exploration globale des données de comptage (par exemple sous forme de carte de chaleur (heatmap)).

Les logiciels d’analyse différentielles (edgeR et DESeq2) ont leur propre normalisation intégrée à la méthode et il n’y a pas besoin d’appliquer de normalisation aux comptages avant un appel à l’une de ces méthodes.

QC (Quality Control)

Comparaison de la taille de librairies

QC (Quality Control)

Comparaison des distributions avant normalisation TMM

QC (Quality Control)

Comparaison des distributions après normalisation TMM

Diagrammes thermiques

  • Appelé également Heatmap
  • Normalisation par ligne ou globale pour l’affichage (Standardisation)
  • Représentation d’une matrice de manière graphique, chaque valeur étant représentée par une cellule colorée.
  • Les colonnes ou les lignes peuvent être ordonnée de plusieurs façons, soit manuellement, soit par une méthode de type clustering.

Diagramme thermique

(Avec standardisation des profils d’expression par z-scores)

Diagramme thermique

(Sans standardisation des profils d’expression par z-scores)

Analyse EdgeR avec statistiques

Diagramme en volcan

  • \(LogFC=1.5, pval=10^{-3}\)

Diagramme en volcan

  • \(LogFC=2.0, pval=10^{-5}\)

Analyse multi-dimensionnelles par Analyse en Composantes Principales

Objectif de l’ACP: Analyse statistique descriptive multi-dimensionnelle

  • Liens entre échantillons: groupes et partitions
  • Liens entre gènes: gènes co-exprimés
  • Liens entre échantillons et gènes
  • Le cas échéant: variables supplémentaires

Structure de base des données

  • En colonnes \(1...K\) variables
  • En lignes \(1...I\) individus

Structure de base de l’exploration par ACP

Représentation d’un nuage de points

Représentation dans \(R^k\). (Figure : F. Husson)

Différentes projections

On cherche à trouver la “meilleure” projection.

C’est une image simplifiée qui “résume” au mieux les données.

Pour cela, on trouve la projection qui restitue au mieux la forme générale du nuage:

(Figure : F. Husson)

Ajustement du nuage des individus

Tableau des températures (F. Husson)

Janv Fevr Mars Avri Mai Juin juil Aout Lati Long Moye Ampl
Bordeaux 5.6 6.6 10.3 12.8 15.8 19.3 20.9 21.0 44.50 -0.34 13.33 15.4
Brest 6.1 5.8 7.8 9.2 11.6 14.4 15.6 16.0 48.24 -4.29 10.77 10.2
Clermont 2.6 3.7 7.5 10.3 13.8 17.3 19.4 19.1 45.47 3.05 10.94 16.8
Grenoble 1.5 3.2 7.7 10.6 14.5 17.8 20.1 19.5 45.10 5.43 10.98 18.6
Lille 2.4 2.9 6.0 8.9 12.4 15.3 17.1 17.1 50.38 3.04 9.73 14.7
Lyon 2.1 3.3 7.7 10.9 14.9 18.5 20.7 20.1 45.45 4.51 11.36 18.6
Marseille 5.5 6.6 10.0 13.0 16.8 20.8 23.3 22.8 43.18 5.24 14.23 17.8
Montpellier 5.6 6.7 9.9 12.8 16.2 20.1 22.7 22.3 43.36 3.53 13.89 17.1
Nantes 5.0 5.3 8.4 10.8 13.9 17.2 18.8 18.6 47.13 -1.33 11.69 13.8
Nice 7.5 8.5 10.8 13.3 16.7 20.1 22.7 22.5 43.42 7.15 14.84 15.2

ACP sur les données Température (Graphe des individus)

ACP sur les données Température (Graphe des variables)

Ajout des variables supplémentaires

Conclusion

L’ACP permet d’obtenir une vision synthétique des données.

  • Groupes d’individus
  • Groupes de variables
  • Quelles variables structurent les données?

-> Passage aux données d’expression RNA-seq.

Exemple sur les données cellulaires précédentes

Cercle des corrélations (gènes)

Analyse fonctionnelle

Après avoir identifié les gènes différentiellement exprimés, nous cherchons à obtenir leur fonction.

Deux manières de procéder:

  • Parcourir la liste et faire une recherche dans la littérature…
  • Utiliser les annotations des gènes pour trouver les fonctions moléculaires ou pathways représentés par ces gènes: -> Faire une Analyse d’enrichissement fonctionnelle.

Analyse fonctionnelle par enrichissement: principe (1)

Analyse fonctionnelle par enrichissement: principe (2)

Analyse fonctionnelle par enrichissement GO

Elles sont basées sur deux composantes:

  • Elles utilisent une Ontologie (vocabulaire contrôlé et stable mis en place par le Gene Ontology Consortium ou autre (Kyoto Encyclopedia of Genes and Genomes, KEGG)).
  • D’où l’utilisation fréquente du raccourcis (Enrichissement GO)
  • Elles sont basées sur un Enrichissement Fonctionnel associé à une validation statistique par Test Hypergéométrique.

Qu’est ce qu’une ontologie ?

Une ontologie est l’ensemble structuré des termes et concepts représentant le sens d’un champ d’informations, que ce soit par les métadonnées d’un espace de noms, ou les éléments d’un domaine de connaissances.

Application au génome: Gene Ontology (Gene Ontology Consortiumhttp://www.amigo.org). 3 ontologies ont été définies.

  • Biological Process
  • Cellular Component
  • Molecular Function

Graphe GO

Exemple d’annotation

RPL35A

Annotations par un vocabulaire contrôlé.

Test d’enrichissement GO

Une catégorie de gènes regroupe n gènes sur le total de N présents sur la puce. La fréquence de départ de cette catégorie est F = n/N.

Ayant obtenus k gènes significativement exprimés ou sous-exprimés, dont p appartiennent à la catégorie C, la fréquence de la catégorie C dans ces gènes est f = k/p.

L’enrichissement est défini comme f/F.

Le test d’enrichissement doit répondre à la question: L’enrichissement est-il statistiquement significatif par rapport à un tirage au hasard? On le fait par Test hypergéométrique.

Exemple de résultat

Rappels sur les extensions de fichiers

  • Fichiers de séquences brutes: .fastq (Compressé en zip: .Fastq.gz)
  • Fichiers de séquences alignées .bam
  • Index de fichiers de séquences alignées .bai
  • Génome complet au format FASTA: .fa
  • Fichier d’expression: Délimité par des tabulations: .txt, .tdf ou .csv.

Rappel des étapes bioinformatiques

  • Contrôle Qualité (FASTQC)

  • Trimming (Trimmomatic)

  • Alignement sur le génome de référence (STAR, aligneur spécifique RNA-seq)

    ou

  • Alignement sur transcriptome de référence(Salmon)

  • Comptage (FeatureCount)

  • Analyse différentielle (EdgeR ou DESeq)

  • Visualisation des données (Morpheus, Bioconductor)

  • Analyse des données (Gene Set Enrichment Analysis)

Rappel sur les différentes normalisations

  • Normalisation sur les comptages (RPKM, TMM, ou autre)
  • Standardisation de l’affichage des heatmaps (affichage centré-réduit par ligne par l’utilisation de z-scores)

Banques de données publiques NCBI

Il s’agit de Dépôts de données liées à des publications répondant à des standards minimum de conservation et de reproducibilité de l’information, et contenant des données brutes de puces à ADN et de NGS.

L’information stockée permet la reproductibilité de l’expérience Leur usage est exigé pour publication (Numéro d’accession).

Dépôts spécialisés:

Bases d’annotations utiles

Autres services NCBI

PubMed:

http://www.ncbi.nlm.nih.gov/pubmed

  • Recherches par auteur, années, titre, contenu du résumé
  • Recherche par publications en lien
  • Possibilité d’ouvrir un compte utilisateur NCBI (Recherches favorites, Recherches automatisées)
  • Un grand nombre de publications sont accessibles de façon gratuite (PubMed Central)

Quelques Programmes utiles

Partie 2: Analyse Single Cell RNA-seq (scRNA-seq)

Partie 2: Analyse Single Cell RNA-seq (scRNA-seq)

  • Etude de l’expression des gènes au niveau de chaque cellule

  • Etudier l’hétérogénéité d’un tissu (Bulk RNA-seq= les expressions sont confondues)

  • Etudier les types cellulaires rares et inconnus

Il est également possible de réaliser d’autres omiques (Atac-Seq, scDNA-seq).

D’après Cao et al. 2017 (Science)

From RNA-seq to Single Cell RNA-seq

Différentes technologies Single Cell

Technologie Nombre de cellules Nombre de paramètres Pour Contre
Cytométrie en flux 1k - 100k 1-15 Simple N Limité
Spectrométrie de Masse 1k-100k 20-50 Plus de paramètres N Limité
RNA FISH ~100 1 Résolution spatiale Complexe, Faible débit
Multiplex FISH 100-aine ~100 Résolution spatiale Complexe à mettre en oeuvre et à analyser
SS2 scRNA-seq 100-1000 ~20000 Haut débit Matrices creuses
Droplets scRNA-seq (Gouttelettes) 100-1e6 ~20000 Haut débit Matrices creuses

Analyse 10X Genomics (Chromium)

  • Le gel contient les billes magnétiques barcodées destinés à identifier chaque cellule de façon unique.
  • Les cellules sont ensuite isolées dans des goutelettes par émulsion grâce à l’apport d’huile. Les cellules associés à chaque goutelette sont barcodées avec le code de la bille correspondante.
  • La technologie 10X Genomics permets la transcription inverse directement dans les goutelettes (mRNA->cDNA) avec l’ajout d’un barcode unique (UMI = Unique Molecular Identifier) par transcrit. Les transcrits (+ barcode + UMI) sont ensuite séquencés.

Schéma Analytique Général

Le travail de bioinformatique se divise en 3 pipelines:

  • Pipeline primaire, intrinsèque au séquenceur: C’est le pipeline qui produit les séquences brutes

  • Pipeline secondaire: Contrôle qualité de base (suppression des doublons), alignement des séquences (10X Genomics: Logiciel CellRanger)

  • Pipeline tertiaire: Analyse fine de l’expression des gènes en fonction de la question biologique (10X Genomics: Logiciel Seurat sous R)

    • Normalisation + Correction Batch (Lots)
    • Estimation de la dimensionnalité des données (PCA)
    • Réduction de dimensionnalité (tSNE, UMAP)
    • Groupement des cellules semblables par clustering
    • Analyse d’expression différentielle

UMIs

  • UMI: Unique Molecular Identifiers (12bp en plus des 16bp des barcodes des cellules = unique par transcrits)
  • Permettent de donner le comptage quasi-exact des molécules dans un séquençage
  • Réduction du bruit d’amplification par déduplication des fragments de lecture (reads)

Détection des cellules de mauvaise qualité

Plot # transcrits vs barcodes: on détecte les cellules dédoublées (barcodes fortement représentés) ou les cellules non viables (barcodes faiblement représentés)

Détection des cellules de mauvaise qualité (2)

Séparation en cellules de haute et basse qualité: Le point d’inflexion et le coude (knee) du graphe précédent représentent la transition entre les deux composants suivants de la distribution totale de comptage:

  • Des gouttelettes vides sans ARN (barcodes représentés par très peu de transcrits, en bas)
  • Des gouttelettes contenant des doublons avec beaucoup d’ARN (barcodes représentés par énormément de transcrits, en haut)

La différence entre ces deux points permets d’avoir une estimation du nombre de cellules dans notre échantillon.

Pipeline Seurat

  • Développé par le labo Satija. (New York University Genome Center)
  • Actuellement disponible en Version 5.0, fonctionne sous R/Bioconductor.
  • Plusieurs tutoriels sont disponibles sur le site des auteurs (https://satijalab.org/seurat/)
  • Il permet de lire les données alignées par le pipeline secondaire CellRanger (10x genomics) et comporte un grand nombre de possibilité d’analyse.
  • Pour la suite, nous nous intéressons au tutoriel basé sur l’étude des cellules circulantes PBMC (Disponible à https://satijalab.org/seurat/v3.2/pbmc3k_tutorial.html )

Seurat: Filtrage basique

  • Filtrage des cellules ayant peu ou pas de doublons
  • Filtrage des cellules ayant de la contamination mitochodriale (Les gènes MT- représentent les gènes mitochondriales)

Filtrage sur les données

Les données ressemblent à
## CD3D  4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2  3 . . . . . 3 4 1 5
## TCL1A . .  . . . . . . 1 . . . . . . . . . . .  . 1 . . . . . . . .
## MS4A1 . 6  . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . .

Elles contiennent principalement des zéros (matrices creuses). On filtre les cellules ayant trop peu de gènes exprimés et les gènes n’étant pas assez exprimés).

Cette caractéristique est utiliseé pour le stockage optimisé de données en mémoire et sur disque.

Normalisation

Méthode utilisée dans le pipeline Seurat:

  • Méthode de mise à l’échelle globale.

  • L’ensemble des valeurs de comptage est divisé par l’expression totale.

  • On multiple le tout par un facteur \(1e4\).

  • Le résultat est normalisé en \(log\).

  • Avant toute réduction de dimensions (ACP ou autre), les données sont centrées reduites (Variance=1, moyenne=0).

Réduction de dimensions linéaire

L’objectif est d’estimer la dimensionalité des données par application d’une ACP (Analyse en composantes Principales) (Ici estimée à 10).

Valeurs propres de la PCA

Les valeurs propres (Eigen Values) de l’ACP permettent l’estimation de la dimensionalité des données.

Clustering

Les cellules sont clusterisées par l’algorithme suivant:

  • Etape 1: Les cellules sont groupées dans un graphe par similarité de distance dans l’ACP par un algorithme KNN (K Nearest Neighbors, méthode des K plus proches voisins). Cet algorithme prend comme paramètre le nombre de clusters à découvrir. On lui donne le nombre de dimensions estimées précédemment par ACP.

  • Etape 2: le clustering obtenu est optimisé par une méthode appelée Algorithme de Louvain.

Réduction de dimensions (UMAP)

Les clusters déterminés à l’étape précédente sont superposés aux plans obtenus par réduction de dimensions UMAP (Uniform Manifold Approximation and Projection for Dimension Reduction)

Analyse d’expression différentielle

L’analyse d’expression différentielle en scRNA-seq se fait par comparaison des groupes (clusters) détectés lors de l’étape de clustering précédente. Elle est identique à ce qui est fait habituellement en RNA-seq classique, et fonctionne par comparaison de l’expression moyenne des cellules entre les deux groupes pondérée par la variance intra-groupe.

Il est donc possible de faire les analyse suivantes:

  • Comparer des clusters deux à deux.
  • Comparer chaque cluster avec le reste des cellules.
  • On peut aussi explorer l’expression des gènes et étudier leur variation entre les clusters.

Exemple de résultat

Comparaison d’expression entre clusters (Violin Plot)

Identification de population

Les populations cellulaires peuvent être identifiées par l’observation de l’expression des marqueurs d’identification.

Et le résultat reporté de façon visuelle sur UMAP

Expression des marqueurs

Il est également possible d’explorer le résultat de la decomposition UMAP en y superposant l’expression de certains gènes.

Etude de plusieurs échantillons

Publication: Steele et al: Multimodal mapping of the tumor and peripheral blood immune landscape in human pancreatic cancer. Nature Cancer 2020.

  • 16 prélèvements d’adénocarcinome pancréatiques.
  • 3 pancreas sains

Analyse simultanée de plusieurs échantillons ?

Intégration de l’ensemble des points dans une matrice unique, normalisation et application d’une réduction de dimensions UMAP.

Décomposition UMAP - Clusters

Décomposition UMAP - Samples

-> Pas de superposition des échantillons - on clusterise les échantillons !

Après intégration par Harmony

Ce programme permet la normalisation globale des données en “recalant” les échantillons.

PDAC vs control comparison

Manual annotation

Dotplot Gene Expression of PDAC samples clusters

Dotplot Gene Expression of Normal samples clusters

Sub-cluster analysis of Macrophages cells

Liens intéressants:

Merci pour votre attention !

Copyright